!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==============================================================================================================================================|
subroutine tvd_mmc (ro, mx, my, mz, bx, by, bz, en, nx, ng, flux_bnd, n_neighs, sys, dx, dt, gm, max_speed)
!==============================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ng, n_neighs, sys
   real(8), intent(in) :: dx, dt, gm, max_speed
   real(8), dimension(n_neighs, sys), intent(inout) :: flux_bnd
   real(8), dimension(nx), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   integer(4) :: i, ip1, im1

   real(8) :: fro(nx), fmx(nx), fmy(nx), fmz(nx), fbx(nx), fby(nx), fbz(nx), fen(nx)
   real(8) :: fro_s(nx, 2), fmx_s(nx, 2), fmy_s(nx, 2), fmz_s(nx, 2), fbx_s(nx, 2), fby_s(nx, 2), fbz_s(nx, 2), fen_s(nx, 2)
   real(8) :: fro_hat(nx), fmx_hat(nx), fmy_hat(nx), fmz_hat(nx), fbx_hat(nx), fby_hat(nx), fbz_hat(nx), fen_hat(nx)
   real(8) :: ro_h(nx), mx_h(nx), my_h(nx), mz_h(nx), bx_h(nx), by_h(nx), bz_h(nx), en_h(nx)
   
   real(8) :: fro_h(nx), fmx_h(nx), fmy_h(nx), fmz_h(nx), fbx_h(nx), fby_h(nx), fbz_h(nx), fen_h(nx)
   real(8) :: fro_hs(nx, 2), fmx_hs(nx, 2), fmy_hs(nx, 2), fmz_hs(nx, 2), fbx_hs(nx, 2), fby_hs(nx, 2), fbz_hs(nx, 2), fen_hs(nx, 2)
   real(8) :: gro(nx, 2), gmx(nx, 2), gmy(nx, 2), gmz(nx, 2), gbx(nx, 2), gby(nx, 2), gbz(nx, 2), gen(nx, 2)
   real(8) :: rro(nx, 2), rmx(nx, 2), rmy(nx, 2), rmz(nx, 2), rbx(nx, 2), rby(nx, 2), rbz(nx, 2), ren(nx, 2)
   real(8) :: hro(nx), hmx(nx), hmy(nx), hmz(nx), hbx(nx), hby(nx), hbz(nx), hen(nx)
   
   real(8) :: rvel, eps

!----------------------------------------------------------------------------------------------------------------------------------------------|
   eps = epsilon (1.0)
   rvel = dt / dx

!----------------------------------------------------------------------------------------------------------------------------------------------|
!  step one
!----------------------------------------------------------------------------------------------------------------------------------------------|      
   call flux (ro, mx, my, mz, bx, by, bz, en, fro, fmx, fmy, fmz, fbx, fby, fbz, fen, nx, gm)
   call flux_split (ro, mx, my, mz, bx, by, bz, en,                                                                                &
      fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                                      &
      fro_s, fmx_s, fmy_s, fmz_s, fbx_s, fby_s, fbz_s, fen_s, nx, max_speed)

   do i = 1, nx-1
      ip1 = i + 1
      fro_hat(i) = fro_s(i, 1) + fro_s(ip1, 2)
      fmx_hat(i) = fmx_s(i, 1) + fmx_s(ip1, 2)
      fmy_hat(i) = fmy_s(i, 1) + fmy_s(ip1, 2)
      fmz_hat(i) = fmz_s(i, 1) + fmz_s(ip1, 2)
      fbx_hat(i) = fbx_s(i, 1) + fbx_s(ip1, 2)
      fby_hat(i) = fby_s(i, 1) + fby_s(ip1, 2)
      fbz_hat(i) = fbz_s(i, 1) + fbz_s(ip1, 2)
      fen_hat(i) = fen_s(i, 1) + fen_s(ip1, 2)
   enddo

   do i = 2, nx
      im1 = i - 1
      ro_h(i) = ro(i) - rvel * (fro_hat(i) - fro_hat(im1))
      mx_h(i) = mx(i) - rvel * (fmx_hat(i) - fmx_hat(im1))
      my_h(i) = my(i) - rvel * (fmy_hat(i) - fmy_hat(im1))
      mz_h(i) = mz(i) - rvel * (fmz_hat(i) - fmz_hat(im1))
      bx_h(i) = bx(i) - rvel * (fbx_hat(i) - fbx_hat(im1))
      by_h(i) = by(i) - rvel * (fby_hat(i) - fby_hat(im1))
      bz_h(i) = bz(i) - rvel * (fbz_hat(i) - fbz_hat(im1))
      en_h(i) = en(i) - rvel * (fen_hat(i) - fen_hat(im1))
   enddo

!----------------------------------------------------------------------------------------------------------------------------------------------|
!  step two
!----------------------------------------------------------------------------------------------------------------------------------------------|
   call flux (ro_h, mx_h, my_h, mz_h, bx_h, by_h, bz_h, en_h, fro_h, fmx_h, fmy_h, fmz_h, fbx_h, fby_h, fbz_h, fen_h, nx, gm)
   call flux_split (ro_h, mx_h, my_h, mz_h, bx_h, by_h, bz_h, en_h,                                                                &
      fro_h, fmx_h, fmy_h, fmz_h, fbx_h, fby_h, fbz_h, fen_h,                                                                      &
      fro_hs, fmx_hs, fmy_hs, fmz_hs, fbx_hs, fby_hs, fbz_hs, fen_hs, nx, max_speed)

   do i = 1, nx - 1
      ip1 = i + 1
      gro(i, 1) = fro_hs(ip1, 1) - fro_s(i, 1)
      gro(i, 2) = fro_s(ip1, 2) - fro_hs(i, 2)
      gmx(i, 1) = fmx_hs(ip1, 1) - fmx_s(i, 1)
      gmx(i, 2) = fmx_s(ip1, 2) - fmx_hs(i, 2)
      gmy(i, 1) = fmy_hs(ip1, 1) - fmy_s(i, 1)
      gmy(i, 2) = fmy_s(ip1, 2) - fmy_hs(i, 2)
      gmz(i, 1) = fmz_hs(ip1, 1) - fmz_s(i, 1)
      gmz(i, 2) = fmz_s(ip1, 2) - fmz_hs(i, 2)
      gbx(i, 1) = fbx_hs(ip1, 1) - fbx_s(i, 1)
      gbx(i, 2) = fbx_s(ip1, 2) - fbx_hs(i, 2)
      gby(i, 1) = fby_hs(ip1, 1) - fby_s(i, 1)
      gby(i, 2) = fby_s(ip1, 2) - fby_hs(i, 2)
      gbz(i, 1) = fbz_hs(ip1, 1) - fbz_s(i, 1)
      gbz(i, 2) = fbz_s(ip1, 2) - fmz_hs(i, 2)
      gen(i, 1) = fen_hs(ip1, 1) - fen_s(i, 1)
      gen(i, 2) = fen_s(ip1, 2) - fen_hs(i, 2)
   enddo

   do i = 2, nx - 1
      im1 = i - 1
      if (abs(gro(i, 1)) .lt. eps) gro(i, 1) = eps
      if (abs(gmx(i, 1)) .lt. eps) gmx(i, 1) = eps
      if (abs(gmy(i, 1)) .lt. eps) gmy(i, 1) = eps
      if (abs(gmz(i, 1)) .lt. eps) gmz(i, 1) = eps
      if (abs(gbx(i, 1)) .lt. eps) gbx(i, 1) = eps
      if (abs(gby(i, 1)) .lt. eps) gby(i, 1) = eps
      if (abs(gbz(i, 1)) .lt. eps) gbz(i, 1) = eps
      if (abs(gen(i, 1)) .lt. eps) gen(i, 1) = eps
      rro(i, 1) = gro(im1, 1) / gro(i, 1)
      rmx(i, 1) = gmx(im1, 1) / gmx(i, 1)
      rmy(i, 1) = gmy(im1, 1) / gmy(i, 1)
      rmz(i, 1) = gmz(im1, 1) / gmz(i, 1)
      rbx(i, 1) = gbx(im1, 1) / gbx(i, 1)
      rby(i, 1) = gby(im1, 1) / gby(i, 1)
      rbz(i, 1) = gbz(im1, 1) / gbz(i, 1)
      ren(i, 1) = gen(im1, 1) / gen(i, 1)
      if (abs(rro(i, 1)) .ge. 1.d0) rro(i, 1) = 1.d0
      if (abs(rmx(i, 1)) .ge. 1.d0) rmx(i, 1) = 1.d0
      if (abs(rmy(i, 1)) .ge. 1.d0) rmy(i, 1) = 1.d0
      if (abs(rmz(i, 1)) .ge. 1.d0) rmz(i, 1) = 1.d0
      if (abs(rbx(i, 1)) .ge. 1.d0) rbx(i, 1) = 1.d0
      if (abs(rby(i, 1)) .ge. 1.d0) rby(i, 1) = 1.d0
      if (abs(rbz(i, 1)) .ge. 1.d0) rbz(i, 1) = 1.d0
      if (abs(ren(i, 1)) .ge. 1.d0) ren(i, 1) = 1.d0

      ip1 = i + 1
      if (abs(gro(i, 2)) .lt. eps) gro(i, 2) = eps
      if (abs(gmx(i, 2)) .lt. eps) gmx(i, 2) = eps
      if (abs(gmy(i, 2)) .lt. eps) gmy(i, 2) = eps
      if (abs(gmz(i, 2)) .lt. eps) gmz(i, 2) = eps
      if (abs(gbx(i, 2)) .lt. eps) gbx(i, 2) = eps
      if (abs(gby(i, 2)) .lt. eps) gby(i, 2) = eps
      if (abs(gbz(i, 2)) .lt. eps) gbz(i, 2) = eps
      if (abs(gen(i, 2)) .lt. eps) gen(i, 2) = eps
      rro(i, 2) = gro(ip1, 2) / gro(i, 2)
      rmx(i, 2) = gmx(ip1, 2) / gmx(i, 2)
      rmy(i, 2) = gmy(ip1, 2) / gmy(i, 2)
      rmz(i, 2) = gmz(ip1, 2) / gmz(i, 2)
      rbx(i, 2) = gbx(ip1, 2) / gbx(i, 2)
      rby(i, 2) = gby(ip1, 2) / gby(i, 2)
      rbz(i, 2) = gbz(ip1, 2) / gbz(i, 2)
      ren(i, 2) = gen(ip1, 2) / gen(i, 2)
      if (abs(rro(i, 2)) .ge. 1.d0) rro(i,2) = 1.d0
      if (abs(rmx(i, 2)) .ge. 1.d0) rmx(i,2) = 1.d0
      if (abs(rmy(i, 2)) .ge. 1.d0) rmy(i,2) = 1.d0
      if (abs(rmz(i, 2)) .ge. 1.d0) rmz(i,2) = 1.d0
      if (abs(rbx(i, 2)) .ge. 1.d0) rbx(i,2) = 1.d0
      if (abs(rby(i, 2)) .ge. 1.d0) rby(i,2) = 1.d0
      if (abs(rbz(i, 2)) .ge. 1.d0) rbz(i,2) = 1.d0
      if (abs(ren(i, 2)) .ge. 1.d0) ren(i,2) = 1.d0
   enddo

   do i = 1, nx
      hro(i) = fro_hat(i) + (rro(i, 1) * gro(i, 1) - rro(i, 2) * gro(i, 2)) * 0.5d0
      hmx(i) = fmx_hat(i) + (rmx(i, 1) * gmx(i, 1) - rmx(i, 2) * gmx(i, 2)) * 0.5d0
      hmy(i) = fmy_hat(i) + (rmy(i, 1) * gmy(i, 1) - rmy(i, 2) * gmy(i, 2)) * 0.5d0
      hmz(i) = fmz_hat(i) + (rmz(i, 1) * gmz(i, 1) - rmz(i, 2) * gmz(i, 2)) * 0.5d0
      hbx(i) = fbx_hat(i) + (rbx(i, 1) * gbx(i, 1) - rbx(i, 2) * gbx(i, 2)) * 0.5d0
      hby(i) = fby_hat(i) + (rby(i, 1) * gby(i, 1) - rby(i, 2) * gby(i, 2)) * 0.5d0
      hbz(i) = fbz_hat(i) + (rbz(i, 1) * gbz(i, 1) - rbz(i, 2) * gbz(i, 2)) * 0.5d0
      hen(i) = fen_hat(i) + (ren(i, 1) * gen(i, 1) - ren(i, 2) * gen(i, 2)) * 0.5d0
   enddo

   do i = 2, nx
      im1 = i - 1
      ro(i) = ro(i) - rvel * (hro(i) - hro(im1))
      mx(i) = mx(i) - rvel * (hmx(i) - hmx(im1))
      my(i) = my(i) - rvel * (hmy(i) - hmy(im1))
      mz(i) = mz(i) - rvel * (hmz(i) - hmz(im1))
      bx(i) = bx(i) - rvel * (hbx(i) - hbx(im1))
      by(i) = by(i) - rvel * (hby(i) - hby(im1))
      bz(i) = bz(i) - rvel * (hbz(i) - hbz(im1))
      en(i) = en(i) - rvel * (hen(i) - hen(im1))   
   enddo

   flux_bnd(1, 1) = hro(ng)
   flux_bnd(2, 1) = hro(nx - ng)
   flux_bnd(1, 2) = hmx(ng)
   flux_bnd(2, 2) = hmx(nx - ng)
   flux_bnd(1, 3) = hmy(ng)
   flux_bnd(2, 3) = hmy(nx - ng)
   flux_bnd(1, 4) = hmz(ng)
   flux_bnd(2, 4) = hmz(nx - ng)
   flux_bnd(1, 5) = hbx(ng)
   flux_bnd(2, 5) = hbx(nx - ng)
   flux_bnd(1, 6) = hby(ng)
   flux_bnd(2, 6) = hby(nx - ng)
   flux_bnd(1, 7) = hbz(ng)
   flux_bnd(2, 7) = hbz(nx - ng)
   flux_bnd(1, 8) = hen(ng)
   flux_bnd(2, 8) = hen(nx - ng)

   return
end subroutine tvd_mmc
